04 - Redução de dimensionalidade e clusterização

Motivação

Em experimentos de transcriptômica de célula única (scRNA-seq), a alta dimensionalidade dos dados — frequentemente composta por milhares de genes — dificulta a análise direta e a interpretação biológica. Para identificar populações celulares de forma robusta, é necessário reduzir essa dimensionalidade preservando as principais variações biológicas. Posteriormente, técnicas de clusterização são empregadas para agrupar células com perfis transcriptômicos semelhantes, revelando estruturas ocultas nos dados, como tipos celulares ou estados funcionais distintos.

Este tutorial descreve as etapas fundamentais de união de amostras, redução de dimensionalidade, clusterização e subclusterização de células, utilizando o pacote Seurat.

Contextualização

Neste tutorial, trabalharemos com objetos do Seurat e utilizaremos as seguintes funções principais:

  • merge() — para unir objetos de diferentes amostras.
  • SCTransform() ou NormalizeData() — para normalização e modelagem da variação técnica.
  • RunPCA() — para a redução inicial de dimensionalidade via Análise de Componentes Principais (PCA).
  • FindNeighbors() e FindClusters() — para identificar agrupamentos celulares.
  • RunUMAP() ou RunTSNE() — para projeção em duas dimensões, facilitando a visualização.
  • Subset() — para subclusterizar grupos de interesse.
library(Seurat)
library(ROGUE)
library(presto)
library(dplyr)
library(here)
library(ggplot2)

knitr::opts_knit$set(
  root.dir = "/home/oandrefonseca/Disciplinas/PPGBM0117.2025.1",
  verbose = FALSE
)

# Aumentar o limite de uso de memoria
options(future.globals.maxSize = 5 * 1024^3) # 5 GB, por exemplo

Unindo amostras

Antes de qualquer redução de dimensionalidade, é necessário unir diferentes amostras em um único objeto Seurat. Para isso, usamos a função merge(). Abaixo um exemplo de codigo de realizar a juncao.

# Listar os arquivos RDS
files <- list.files(
  path = "path/to/your/objects", pattern = "\\.rds$", full.names = TRUE)

# Carregar todos os objetos em uma lista
seurat_list <- lapply(files, readRDS)

# Criar identificadores a partir dos nomes dos arquivos
sample_ids <- gsub("\\.rds$", "", basename(files))

# Fazer o merge
seurat_merged <- merge(
  x = seurat_list[[1]],
  y = seurat_list[-1],
  add.cell.ids = sample_ids,
  project = "ProjectName"
)

Salvando o objeto do Seurat

if(!file.exists("./data/seurat_object.RDS")) {
  
  saveRDS(file = "./data/seurat_object.RDS", object = seurat_merged)
  
}
# Objeto contendo todas as amostras
seurat_object <- readRDS(file = here::here("./data/seurat_object.RDS"))

Normalização

seurat_object <- NormalizeData(
  seurat_object, normalization.method = "LogNormalize", scale.factor = 10000)
Normalizing layer: counts
seurat_object <- FindVariableFeatures(
  seurat_object, selection.method = "vst", nfeatures = 2000)
Finding variable features for layer counts
seurat_object <- ScaleData(
  seurat_object, features = VariableFeatures(seurat_object))
Centering and scaling data matrix
Warning: Different features in new layer data than already exists for
scale.data

Redução de dimensionalidade (PCA)

seurat_object <- RunPCA(seurat_object, npcs = 50, verbose = FALSE)
ElbowPlot(seurat_object, reduction = "pca")
Figure 1
Note

O ElbowPlot(seurat_object) é utilizado para identificar o número apropriado de componentes principais (PCs) a serem considerados nas análises downstream, como clusterização e redução de dimensionalidade. O gráfico exibe a variação explicada por cada PC; o ponto onde há uma inflexão (“cotovelo”) indica que PCs adicionais contribuem pouco para a variabilidade dos dados. Selecionar o número correto de PCs é essencial para capturar a estrutura biológica relevante, minimizando ruído técnico.

Visualizando dimensões

# Examine and visualize PCA results a few different ways
VizDimLoadings(seurat_object, dims = 1:2, reduction = "pca")
Figure 2

Extraindo informacões do PCA

# Get cell embeddings and feature loadings stored on pbmc[['pca']]@cell.embeddings
Embeddings(seurat_object, reduction = "pca")

# stored in pbmc[['pca]]@feature.loadings
Loadings(seurat_object, reduction = "pca")

Construção da matriz de vizinhança

seurat_object <- FindNeighbors(seurat_object, dims = 1:20)
Computing nearest neighbor graph
Computing SNN
Note

A função FindNeighbors constrói uma matriz de vizinhança baseada nos embeddings de redução de dimensionalidade (tipicamente PCA). Cada célula é conectada às suas células mais próximas no espaço de baixa dimensão, formando um grafo de similaridade.

Esse grafo é o alicerce para a clusterização subsequente (FindClusters()), permitindo identificar comunidades de células com perfis transcriptômicos semelhantes. O parâmetro dims especifica quais componentes principais serão usados para calcular as distâncias entre as células.

Clusterização

seurat_object <- FindClusters(
  seurat_object, algorithm = 4, resolution = c(0.25, 0.5, 1.0))
Note

A função FindClusters realiza a clusterização das células com base na matriz de vizinhança previamente calculada.
Ela aplica algoritmos de otimização de comunidades (como o algoritmo Louvain ou Leiden) para agrupar células com perfis transcriptômicos semelhantes.

O parâmetro resolution controla a granularidade da clusterização, por exemplo: - Valores mais baixos (< 0.5) resultam em menos clusters (agrupamentos maiores). - Valores mais altos (> 0.5) resultam em mais clusters (agrupamentos menores e mais específicos).

A escolha do resolution deve ser guiada tanto por critérios biológicos quanto por inspeção visual dos agrupamentos gerados.

# Explore as resolutions
DT::datatable(seurat_object@meta.data)
Warning in instance$preRenderHook(instance): It seems your data is too big for
client-side DataTables. You may consider server-side processing:
https://rstudio.github.io/DT/server.html

Projeção com UMAP

seurat_object <- RunUMAP(seurat_object, dims = 1:20)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
19:29:32 UMAP embedding parameters a = 0.9922 b = 1.112
19:29:32 Read 67633 rows and found 20 numeric columns
19:29:32 Using Annoy for neighbor search, n_neighbors = 30
19:29:32 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
19:29:39 Writing NN index file to temp file /tmp/Rtmpl4WPQu/file1701329171d72
19:29:39 Searching Annoy index using 1 thread, search_k = 3000
19:30:03 Annoy recall = 100%
19:30:04 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
19:30:05 Initializing from normalized Laplacian + noise (using RSpectra)
19:30:23 Commencing optimization for 200 epochs, with 3039738 positive edges
19:30:48 Optimization finished
# Por padrao apenas os clusters da ultima resolucao sao apresentados
DimPlot(seurat_object, reduction = "umap", label = TRUE)
Figure 3
Note

Importante: A clusterização de scRNA-seq não é um processo exato. Ela deve ser encarada como uma ferramenta de geração de hipóteses. Sempre valide biologicamente seus agrupamentos antes de tirar conclusões definitivas.

Expressão diferencial

# O parametro `max.cells.per.ident` foi acrescentado para otimizar o tempo computacional. Remove-o para analises reais.
cluster_markers <- FindAllMarkers(seurat_object, max.cells.per.ident = 1000)
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4
Calculating cluster 5
Calculating cluster 6
Calculating cluster 7
Calculating cluster 8
Calculating cluster 9
Calculating cluster 10
Calculating cluster 11
Calculating cluster 12
Calculating cluster 13
Calculating cluster 14
Calculating cluster 15
Calculating cluster 16
Calculating cluster 17
Calculating cluster 18
Calculating cluster 19
Calculating cluster 20
Calculating cluster 21
Calculating cluster 22
Calculating cluster 23
Calculating cluster 24
Calculating cluster 25
Calculating cluster 26
Calculating cluster 27

Inspecionando DEGs

# Removendo resultados nao significativos
cluster_markers_top <- cluster_markers %>%
  filter(p_val_adj <= 0.01) %>%
  group_by(cluster) %>%
  top_n(20, avg_log2FC)

DT::datatable(cluster_markers_top)

Salvando tabela dos DEGs

saveRDS(cluster_markers, file = "./data/project_differential.RDS")

Analise de pureza dos clusters

downsampled_object <- subset(x = seurat_object, downsample = 1000)
counts_matrix <- LayerData(
  downsampled_object, assay = "RNA", layer = "counts")

Calculando entropia

entropy_calculation <- SE_fun(counts_matrix)
rogue_value <- CalculateRogue(
  entropy_calculation, 
  platform = "UMI")

rogue_value
[1] 0.2160008
rogue_result <- rogue(
  counts_matrix, 
  labels = downsampled_object@meta.data$seurat_clusters, 
  samples = downsampled_object@meta.data$orig.ident, 
  platform = "UMI", 
  span = 0.6
  )
rogue.boxplot(rogue_result)
Figure 4